Draft version February 8, 2012 

Preprint typeset using I^T^X style cmulatcapj v. 5/2/11 



o 

IX, 



< 

6 

(N 
> 

in 

o 
in 

en 

r^ 

o 
o 



X 



FEEDBACK FROM CENTRAL BLACK HOLES IN ELLIPTICAL GALAXIES: TWO-DIMENSIONAL MODELS 

COMPARED TO ONE-DIMENSIONAL MODELS 

Gregory S. Novak, 1 Jeremiah P. Ostriker, 1 , and Luca Ciotti 2 

Draft version February 8, 2012 

ABSTRACT 

We extend the black hole (BH) feedback models of Ciotti, Ostriker, and Proga to two dimensions. In 
this paper, we focus on identifying the differences between the one-dimensional and two-dimensional 
hydrodynamical simulations. We examine a normal, isolated L* galaxy subject to the cooling flow 
instability of gas in the inner regions. Allowance is made for subsequent star formation, Type la and 
Type II supernovae, radiation pressure, and inflow to the central BH from mildly rotating galactic gas 
which is being replenished as a normal consequence of stellar evolution. The central BH accretes some 
of the infalling gas and expels a conical wind with mass, momentum, and energy flux derived from both 
observational and theoretical studies. The galaxy is assumed to have low specific angular momentum 
in analogy with the existing one-dimensional case in order to isolate the effect of dimensionality. The 
code then tracks the interaction of the outflowing radiation and winds with the galactic gas and their 
effects on regulating the accretion. After matching physical modeling to the extent possible between 
the one-dimensional and two-dimensional treatments, we find essentially similar results in terms of 
BH growth and duty cycle (fraction of the time above a given fraction of the Eddington luminosity). 
In the two-dimensional calculations, the cool shells forming at 0.1-1 kpc from the center are Rayleigh- 
Taylor unstable to fragmentation, leading to a somewhat higher accretion rate, less effective feedback, 
and a more irregular pattern of bursting compared to the one-dimensional case. 



1. INTRODUCTION 

Nearly all massive galaxies are thought to harbor su- 
permassive black holes (SMBHs), and the properties 
of the black holes (BHs) are known t o be well corre- 
lated with those of their h ost galaxies (Gcbhardt et all 
2000: Ferraresc & Mcrritt 2000; Tremaine et al. 2002; 
Novak et alj2006tlGultekin e t al. 200jJ. The direction of 



the causal link between SMBHs and host galaxy proper- 
ties remains unclear, but the existence of such a correla- 
tion suggests an intimate connection between the physics 
of BH growth and galaxy formation. There has been in- 
tense theoretical interest in the ability of active galactic 
nuclei (AGNs) to affect the star formation histories and 
observed colors of galaxies (e.g.. iDi Matteo et al.l 120051 ; 
iCroton et al.ll2006h . There is certainly sufficient energy 
and momentum available to have a profound effect on 
the structure and dynamics of the interstellar medium 
(ISM). 

Many physical processes can act to couple accretion by 
the central BH to the galaxy's ISM. The most obvious 
and well-observed result of SMBH accretion is the prodi- 
gious luminous output. Radiatively efficient accretion 
converts a significant fraction of the rest mass energy of 
the infalling gas to radiation, and the bulk of the matter 
that enters BHs over the history of the uni verse is know n 
to do so in a radiatively efficient manner (|Sortanlll982f ). 
The emitted photons impart energy and momentum to 
the galaxy's ISM via electron scattering, photoionization, 
scattering due to atomic resonance lines, and absorption 
by dust grains. 

Accreting SMBHs are also known to drive broad ab- 
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sorption line (BAL) winds that convey mass, momen- 
tum, and energy to the surrounding galaxy. These 
winds are launched by electromagnetic processes within 
a few hundred gra vitational radii of the central BH (e.g., 
iProga et alll2000l ). but by the time the wind galactic 
length scales, the energy is carried by the kinetic motion 
of the gas. Only fraction of the energy output of the 
AGN will be converted to a wind, that wind will pro- 
ceed to inject nearly all of its energy and momentum into 
the galaxy's ISM. This high interaction efficiency has in- 
spired strong theoretical interest in the effect of mechani- 
cal AGN feedback on galaxy propertie s (Di Matte o et al.l 
[20051 : ISprmgel et al J 120051: 1 Johansson et alJl2009f )~~ 

The combined effect of these feedback processes can 
dramatically affect the galaxy and act to reduce subse- 
quent SMBH accretion. Physical modeling of the com- 
plex processes involved is beginning to reach a level where 
observational tests of the models are possible. Building 
on pr evious work in one dimension (ICiotti fc Ostriker 
" ICiotti et al.l 120091* IShin et all 1201ft 
we perform two-dimensional simula- 
tions of SMBH accretion in the context of normal but 
isolated elliptical galaxies. Thus, the source of accreted 
gas is the fraction (~ M*/6) ejected into the ISM due 
to normal processed of stellar evolution, rather than the 
comparable amounts of gas added via cosmological in- 
fall. We take care to resolve the inner length scales 
of accretion (the Bondi radius) while also running the 
simulations long enough to resolve galactic length and 
timescales. Many physical processes must be included 
and numerically a large dynamic range is needed. Al- 
though quite well developed from the point of view of 
the included physics, our previous one-dimensional sim- 
ulations obviously could not address the impact of gen- 
uinely two-dimensional phenomena on the feedback phe- 
nomenon. We recall here the most important of them, 
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that will be discussed in this paper. 

The first one concerns the possible instabilities (e.g., 
Rayleigh-Taylor, RT) that can affect the evolution of 
the cold shells that appear during the evolution of one- 
dimensional models. The central bursts are caused by 
these shells, and the shells have a two-fold origin. At the 
beginning o f each majo r burst, the classical Field cooling 
instability (|Fieldl [l965t ) appears around 1 kpc from the 
center, due to the local critical balance between heating 
and cooling. As the density increases, the shell starts to 
fall toward the center and compresses the gas. When a 
burst first appears, shock waves are sent from the cen- 
ter toward the falling shell, and a series of sub-bursts 
and consequent reflected shock waves impact on the cold 
shell, increasing its density still further. As the bulk of 
star formation in the one-dimensional models happens 
in these cold shells, it is very important to understand 
the cold shell physics, not only from the point of view 
of central accretion, but also for the starburst which oc- 
curs within the cold shell. From the short description 
above, a few obvious questions arise: for example, will 
the allowance of the additional degree of freedom still 
lead to a formation of a cold shell near the center in the 
case of an aspherical galaxy? What is the effect of non- 
zero angular momentum in the gas? Will the cold shell 
fall toward the center as in the one-dimensional simula- 
tions or it will break up due to (RT) instability? Even 
more important, what is the fate of the multiple inter- 
acting shocks? Will the accretion still be characterized 
by strong bursts separated by long time intervals or will 
the breakup of the shells lead to cold fingers of dense gas 
being accreted in a more or less steady flow while hot gas 
flows outward, therefore resulting in flows that at each 
radius are partially accreting and partially outflowing? 

The second reason to move to two-dimensional simu- 
lations is to explore the interaction (and the consequent 
mechanical feedback) of the conical nuclear wind with 
the galaxy ISM. In our previous one-dimensional sim- 
ulations this interaction was necessarily described as a 
spherical average of an inherently non-spherical effect, 
even though we had taken into account several physical 
aspects of the phenomenon via a time-dependent differ- 
ential equation. Clearly, a two-dimensional simulation 
is also needed to explore Kelvin-Helmholtz instabilities 
at the interface between the outflowing conical wind and 
the ISM. 

In the present work, we focus on two-dimensional sim- 
ulations of a galaxy with very low specific angular mo- 
mentum in order to make close contact with the existing 
one-dimensional simulations. We would like to isolate 
the effect of increasing the dimensionality of the simu- 
lation. The adopted angular momentum profile is con- 
sistent with the slowe st of the SAURON slow- rotators 
(jEmsellem et al.ll2004| ). In future work, we will expand 
our treatment to include angular momentum transport 
via t he standard a prescription (|Shakura fc Sunvaevl 
119731) a nd more recent models bas ed on gravitational 
torques (|Hopkins &: Quataertll2010al ). 

There have been many numerical simulations of SMBH 
accretion and the subsequent effects on the galaxies 
containing the resulting AGN. Nearly all of the efforts 
to date can be classified into three broad categ ories. 
Pi Matteo et al . (2005). lDebuhr et al.l (|201d . 120111) . and 
Johansson et al.l (|2009Q are examples where the simu- 



lations cover length scales from ~ 100 pc to tens of 
kpc and timescales from a fraction of a Myr to several 
Gyr. Galactic length and timescales are resolved, but 
the SMBH accretion and feedback processes are con- 
sid ered to be sub-reso l ution. Complementary studies 
bv IKurosawa &: Progal (|2009af ) and iKurosawa fc Progal 
(2009rJ are examples of multi-dimensional simulations 
that cover the length scales from a few AU to ~ 1 pc. 
Length and timescales relevant to SMBH accretion are 
resolved, and the generation of radiatively driven winds is 
computed, but these simulations do not approach galac- 
tic leng th or timescales, and infall rates are taken as 
given iHopkins fc Quataertl (|2010bD and iLevine et al.l 
(|2008f) are examples of a multi-resolution studies of 
SMBH accretion involving progressively higher spa- 
tial resolution simulations run for progressively shorter 
times. The highest spatial resolution simulations go 
down to a fraction of a pc and are run for about one 
Myr of simulation time. These simulations spatially re- 
solve the accretion process, but do not reach galactic 
timescales. Therefore, they cannot self-consistently cal- 
culate the effect of AGN feedback on the gas in the galaxy 
as a whole and the subsequent SMBH accretion. 

Finally, there have been several numerical studies 
of accretion by intermediate-mass black holes (IMBHs) 
with the goa l of understanding BH growth in th e 
early universe (jAlvarez et al.l200a|Park fc Ricottil2010l) . 
In terms of dimensionless length scales r/rBondi or 
?V r 'Schwarzschiid, some of these simulations are similar to 
our simulations. However, studies of IMBH accretion 
focus on BHs with very small masses (100-1000 M Q ) 
compared to those presented here. Therefore, the rele- 
vant physical length and timescales are much smaller and 
the physical sources of the infalling gas are very different 
from those considered in the present work. 

Our goal is to resolve both the relevant accretion 
length and timescales while at the same time resolv- 
ing galactic length a nd timescales (cf. ILevine et al.ll2008t 
lAlvarez et al.l 120091 ) . There have been only a few at- 
tempts to perform multi-dimensional simulations that 
bridge the gap between galactic and SMBH scales, 
although several papers have examined the interac- 
tion between an outflowing wind/jet with specified 
prop erties and the su rrounding intcrgalactic medium 
(cf. iMetzler fc Evrardl 11994 lOmma fc Binnevl l2004t 
i Siiacki et al.ll2007t [Sternberg fc Sokerl2008HReeves et all 
2009; Fa bian et al.ll2009HArieli et aLlboidf ) 

The present work is an attempt to simultaneously re- 
solve the inner length scales relevant to SMBH accretion 
(a few pc), outer length scales relevant to galaxies (tens 
of kpc), inner timescales relevant to SMBH accretion (a 
few years) , and outer timescales relevant to galaxies and 
stellar evolution (10 Gyr). However, the region inside 
of 1 pc including the disk and the SMBH itself are still 
treated as sub-resolution physics and we compute the 
output from these regions as time-dependent functions 
of the input to them, utilizing formulae from the above 
quoted sources. 

We take particular care to resolve the inner scales 
where the rate of accretion is set (the Bondi radius) even 
for the hot gas in the system. This is important because 
heating due to AGN radiation or mechanical energy in- 
put may raise the gas temperature as it falls toward the 
SMBH, decreasing the relevant Bondi radius. For exam- 
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pie, lAlvarez et al.l (|2009f) performed cosmological simu- 
lations of accretion onto the SMBH left by the first stars 
at high rcdshift. Their resolution was sufficient to resolve 
the Bondi radius for the cold gas in their simulation, but 
not for the much hotter gas affected by feedback. This 
means that if the heating would have been effective be- 
tween the resolution of their simulation and the actual 
hot-gas Bondi radius, then in reality the gas would never 
have made it to the SMBH. It is important to realize 
that the material emitted by the evolving stars (the fuel 
for SMBH accretion) is thermalized by stellar velocity 
dispersion to the equivalent temperature, and that the 
stellar velocity dispersion near the BH increases (as 1/r 
in the isotropic case). Therefore, the gas is injected at 
higher and higher temperature approaching the central 
BH and the location of the Bondi radius does not change 
significantly even if the SMBH grows in mass. Care must 
be exercised in the numerical treatment of the central re- 
gions. To put it simply, Bondi-type accretion cannot be 
computed properly unless the detailed thermal state of 
the gas is followed to well within the radius at which 
gravitational and thermal energies are balanced. 

Section [2] describes the simulations and our basic ap- 
proach, Section [3] gives our results, and Section 0] sum- 
marizes our conclusions. 

2. METHODS 

We use th e Zeus hydrodynamics code 

IjStone fe Norman] I1992D because of its ability to 
utilize spherical geometry and variable cell sizes. The 
code is not adaptive in the sense that it does not 
dynamically decide where to place extra cells, but 
the variable cell sizes allow us to identify a center of 
interest and take advantage of increased resolution near 
that point. We extend the code with mass, energy, 
and momentum source terms appropriate for stellar 
evolution in elliptical galaxies, SMBH accretion, and 
feedback resulting fro m the accr e tion, f ollowing the 
treatment described in lCiotti et al.l (|2009bl ). 

The simulation grid is the meridional plane in spher- 
ical coordinates where all quantities are assumed to be 
axisymmetric. The radial bins are the same as in the 
one-dimensional simulations: 120 cells covering 2.5 pc to 
250 kpc where each cell is 10% larger than the previous 
cell. Requiring the cells to have an aspect ratio of unity 
gives 30 angular cells. We exclude the region from to 
0.05 radians near either pole because cell volumes go to 
zero due to the coordinate singularity there. 

The simulations are not limited by the total volume 
of computation required, but rather by the time to solu- 
tion. The Courant-Friedrich-Levy (CFL) condition for 
the central grid cells is quite severe: time steps are ~ 
10 years, and we would like to run the simulation for 12 
Gyr. Numerical schemes which allow the time step to in- 
crease with increasing radius could be implemented, but 
the estimated speedup is quite limited. For a single pro- 
cessor this would save a factor of N r (l~(r i /r ) 1 / Nr ) ~ 10 
where N r is the number of grid cells in the r-direction, 
Ti is the innermost radius, r is the outermost radius, we 
have assumed log-spaced bins in radius, and the numeri- 
cal value is for our chosen parameters otri/r — 10~ 5 and 
N r = 120. In a multi-processor environment, the reduc- 
tion in wallclock-time is given by the same formula with 
ri,r , and N r referring to the cells on a given processor 



rather than the whole simulation. For our simulations 
running on eight processors, the speedup would be only 
-40%. 

As the mass of the central SMBH increases, one might 
think that the Bondi radius increases so that a moving 
inner grid could help in speeding up the simulations. As 
discussed above, this will not work since the gas losses 
from the evolving stars near the SMBH are thermalized 
at the local value of the velocity dispersion, which in- 
creases with BH mass. In practice, the location of the 
"Bondi radius" remains approximately fixed for the en- 
tire time of the simulation, and varies chiefly due to fluc- 
tuations in the local speed of sound caused by feedback 
from the SMBH and by the accretion of cold gas. 

Our requirements for the spatial resolution near the 
center of the simulation, the amount of time for which 
the simulation must run, and the CFL condition near the 
center have constrained us to use a rather small number 
of cells. We have performed a few simulations where the 
number of grid cells is doubled in each dimension, and the 
simulation ran for ~ 1.3 Gyr rather than 12 Gyr. Most 
physical quantities did not change significantly. The ex- 
ception was the change in the BH mass, which was a 
factor of two larger. Smaller cells allow denser blobs and 
filaments to form, which have an easier time making their 
way to the center of the simulation. It is also important 
to note that the size of the cells grows linearly with radius 
so that physical structures of constant size (e.g., molec- 
ular clouds) will be resolved progressively more poorly 
further from the center of the simulation. 

2.1. Physics 

For the most part the input phy sics is the same as 
that described in iCiotti et al.l ([2010) with a few excep- 
tions described in detail belo w. A complete descriptio n 
of the input physics is gi ven inlCiotti fc Ostrikerl (|2007t l. 
ICiotti et alJ (J2009bQ . and lSazonov et al.l (|2005ir Here we 
repeat the most important aspects. 

The total gravitational potential of the model galaxy is 
assumed to be a singular isothermal sphere plus a point 
mass for the central BH. This is good agreement with re- 
cent obse rvations of the total mass profile of early-type 
galaxies ()Gavazzi et al.l 120071 I2008D . For simplicity we 
maintain this model in to the smallest radii, although 
more complicated models may be more appropriate in- 
side of a fraction of the half-light radius. The velocity 
dispersion parameter of the isothermal potential is 260 
km s" 1 . The gas is not self-gravitating. 

The stellar distribution is given by a Jaffe profile with 
a total mass of 3 x 10 u Mq and a projected half-mass 
radius of 6.9 kpc. The mass-to-light ratio is assumed to 
be spatially constant and is equal to 5.8 in solar units in 
the B band at the present time. 

Energy and mass input due to stellar evolution, type 
la supernovae, star formation, type II supernovae are 
exactly the same. The specific energy of the material 
provided by stellar evolution is calculated according to 
the solution of the Jeans equation for the given stellar 
de nsity distribution an d gravitational potential as given 
in lCiotti et al.l (|2009aH . Gas heating and cooling due to 
atomic processes as well as Compton heating and bound- 
free absorption of radiation f rom the AGN are ca lcu- 
la ted using the form ulae from iSazonov et al.1 ([2005D , as 
in lCiotti et al.l (|2010f ). 
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Gas temperatures are bounded from below by the fact 
that our atomic cooling curve has an exponential cutoff 
below 10 4 K and from above by our assumed broad-line 
wind velocity of 10,000 km s _1 , corresponding to 2.5 xlO 9 
K if all of the kinetic energy is converted to thermal en- 
ergy. However, this very hot gas will always be strongly 
outflowing until it distributes its energy to a much larger 
mass of gas so that the temperature is of order the virial 
temperature of the halo. The upper temperature bound 
for infalling gas is given by the maximum temperature to 
which the AGN photons can heat the gas. This is given 
by the Compton te mperature, which w e take to be 2.1 
keV = 2.5 x 10 7 K (|Sazonov et alJl2005h . 

Star formation is implemented using the standard 
Schmidt-Kennicut prescription. The star formation rate 
density is given by 
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where p is the local gas density, 7?form = 0.1 is a dimen- 
sionless parameter for the rapidity of star formation and 

(2) 
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where E is the internal energy per unit volume of the 
gas, C is the volumetric cooling rate, G is the Newto- 
nian gravitational constant, r is the distance from the 
center of the galaxy, and v c (r) is the circular velocity 
as a function of radius. Type II supernovae then return 
mass and energy to the ISM such that the mass returned 
is 20% of that processed through stars and the energy is 
Esn = 4 x 10~ 6 mgFC 2 . These values for the energy and 
mass returned result from our assumption of a Salpctcr 
initial mass function with a low- mass cutoff of 0.1 Mq, 
together with our assumption that every star above 8Mq 
injects its entire mass and 10 51 erg into the ISM. 

In the momentum equation, the force exerted on gas 
due to Compton scattering of photons from the AGN 
is treated as in the one-dimensional code. Note that 
because the vast majority of the photons from the 
AGN have energies far below the electron rest mass en- 
ergy, Compton scattering takes place in the coherent 
limit where photon energies and numbers are preserved. 
Therefore, in both the optically thin and optically thick 
limits, the force per unit mass on a fluid clement is 



F/m 






(7) 



where «es ls the electron scattering opacity and L(r) is 
the total luminosity emitted interior to r (that is, there 
is no e~ T factor, where t is the optical depth). 

The force exerted by absorption of AGN photons by 
atomic lines is also the same as in the one-dimensional 
code: dp/ dt — l/c, where I is the energy absorbed per 
unit time in a given cell. In the work reported on in 



this paper, the radiation momentum associated with dust 
absorption (cf. iDebuhr et al.ll20lTI ) is not included in the 
two-dimensional calculation although it is included in the 
one-dimensional treatment. 

Some physical processes are modeled using the same 
basic assumptions as in the one-dimensional code, but 
two-dimensional simulations allow us to improve the im- 
plementation. 

This is the case with the BAL wind. The one- 
dimensional and two-dimensional simulations use the 
same formulae to compute the total mass, momentum, 
and energy injected into the ISM for a given SMBH ac- 
cretion rate. However, the one-dimensional simulations 
use a prescription based on pressure balance between the 
outgoing wind and the ambient gas in the galaxy in or- 
der to compute how the mass, energy, and momentum 
are radially distributed. In the two-dimensional simula- 
tions, we can adopt the more realistic treatment of simply 
injecting the desired mass, energy, and momentum into 
the innermost radial cells and self-consistently comput- 
ing the radial transport of these quantities. 

The one-dimensional code specifies the opening angle 
of the wind for the purpose of solving for the radial dis- 
tribution of the deposited material. The two-dimensional 
code requires a more precise specification of the depen- 
dence of the outflowing material on angle from the pole. 
We choose the simple parameterization: 



dq _ (w + 1) sin6>|cos"6>|Q 

~dl ~ 4~t 



(8) 



where q is a conserved quantity (mass, energy, or ra- 
dial momentum) , Q is the total amount of the conserved 
quantity to be injected. We solve for n such that the an- 
gle 9 enclosing half of the conserved quantity is the same 
as the opening angle specified in the one-dimensional 
code. For the fixed-opening-angle models considered 
here, n = 2, so that the half-opening angle enclosing 
half of the energy is ~ 45°. In terms of solid angle, this 
means that the wind is visible from ~ 1/4 of the avail- 
able viewing angles. This fraction is in agreement with 
observations of the fraction of obscured and unobscured 
AGNs under the assumption that the two populations 
are made up of a single population of objects that differ 
only in viewing angle. 

Finally, some processes of minor importance have been 
omitted from the two-dimensional simulations so far. 
This includes the prescription for star formation and stel- 
lar remnant formation in the central sub-grid SMBH ac- 
cretion disk. The two-dimensional code implements no 
star formation in the disk — all gas that flows in eventu- 
ally either flows into the SMBH or back onto the simu- 
lation grid as a conical wind. 

The major difference between the calculations pre- 
sented here and the one-dimensional models is that the 
two-dimensional models do not yet implement any op- 
tically thick radiative transfer algorithms. We assume 
that the gas is always optically thin at all wavelengths at 
all times. There is also no treatment of the temperature 
dependence of dust opacity. Finally, the radiation pres- 
sure due to stars formed in the simulation is also omitted 
from the momentum equation. 

Experiments using one-dimensional models have shown 
that omitting the radiation pressure on dust does not 
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greatly affect the overall SMBH mass growth when these 
other feedback mechanisms are included (less than a fac- 
tor of two). However, star formation can be affected 
because infalling cold shells can have their falling time 
extended by radiation pressure on dust, allowing more 
stars to form as the gas is falling. Improving the treat- 
ment of radiative transfer in the two-dimensional code 
will be the subject of future work. 

Computational efficiency requires that we limit the ve- 
locities and sound speeds in the simulation to prevent the 
CFL condition from restricting the time step too dramat- 
ically. At the high end, we impose a limit of 30,000 km 
s~ 4 on the gas velocities and sound speeds. The source 
term with the highest specific energy is the mechanical 
wind due to SMBH accretion, which injects gas at 10,000 
km s^ 1 . The limit is three times higher, corresponding 
to 10 times the energy density, so this limit should not 
greatly affect overall evolution of the system. 

In the one-dimensional simulations, we imposed the 
commonly used requirement that the ISM temperature 
be above 10 4 K. This was necessary because the one- 
dimensional code uses a fully explicit time integration 
scheme and low thermal energies imply exceedingly short 
cooling times. In the two-dimensional calculations, we 
use a semi-implicit time integration method, making it 
possible to set the lower bound on temperature at a lower 
value. We impose three limits dictating that the gas does 
not drop below the temperature of the cosmic microwave 
background, the effective temperature associated with 
the AGN radiation field or the effective temperature of 
the stellar radiation field. The two-dimensional simu- 
lations make use of the same cooling curve as the one- 
dimensional simulations, which has an exponential cutoff 
below 10 4 K. Therefore, gas cannot reach these low tem- 
peratures by cooling alone, but must undergo a dramatic 
expansion. This occasionally occurs due to the violent 
gas motions near the center of the simulation. 

We calculate the temperature of the cosmic microwave 
background assuming an Einstein-de Sitter universe with 
an age of 13.7 Gyr where the beginning of the simulation 
corresponds to 1.7 Gyr after the big bang. The difference 
between an Einstein-de Sitter universe and the standard 
ACDM cosmology is minimal for our present purpose. 

The effective temperature of the photons from the 
AGN is just 

( L V /4 



T = 



where a = 8tt 5 k B / 15c 3 h 3 is the radiation constant. 

It is easy to show that for a spherically symmetric dis- 
tributed luminosity source, the energy density in the pho- 
ton field at R is 



e{R) 



j(r) dr , R 



where j(r) is the luminosity per unit volume. 
However, we adopt the expression 



e(R) = 



2 j(r) dr 
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because it gives simple analytic integrals and has the cor- 
rect form in the limits when r S> R and r <C R. The inte- 
grand of Equation (fTTj) underestimates the integrand of 



Equation (|10[) when r is near R. However, the difference 
between the two expressions is less than 30% provided 
that r and R differ by at least a factor of two. Given 
that l(r) is typically a steeply falling function of r, the 
integral is not typically dominated by the region where 
r~R. 

For the present case of a Jaffe profile assuming a con- 
stant stellar mass-to-light ratio, this gives 
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where y = R/r*, r* is the scale radius from the Jaffe 
profile, L is the total luminosity, and the expression has 
a finite limit when y = 1. 

There are also some minor computational differences 
between the two codes. The one-dimensional simulations 
use a fully explicit algorithm for time integration. In the 
two-dimensional code, the time integration is mostly ex- 
plicit, but we use a semi- implicit scheme for radiative 
cooling because the cooling timescales can become very 
short. If the time step required by the CFL condition 
is shorter than the cooling time, the code updates gas 
temperatures due to radiative cooling explicitly. If the 
cooling time is shorter than the CFL time step, the code 
evaluates the cooling function at the advanced time, us- 
ing a bisection algorithm to find the root of the resulting 
equation for the gas temperature at the advanced time. 

2.2. AGN Feedback Parameterization 

The BAL originating near the SMBH provides energy, 
momentum, and mass from a wind to the ISM according 
to the equations: 



Mbh = MfafcuAl + V) , 


(13) 


n = 2e w c 2 /v 2 v , 


(14) 


% = cwMbhc , 


(15) 


Pw = 2e w c 2 M B K/vw , 


(16) 



and 



M w = 2e w c 2 M BK /v 



W, (17) 

where vw is the velocity of the BAL, ta ken here to be 
10,000 km s _1 . As discussed in detail in lOstriker et al.l 
(2010), these expressions guarantee that the mass, en- 
ergy, and momentum carried by the wind are self- 
c onsistent. 

iCiotti et al.l ()2009bl ) considered two classes of mechan- 
ical feedback, denoted the A and B models. For the A 
models, both the mechanical efficiencies and the wind 
opening angles were independent of the accretion rate. 
The B models have mechanical efficiencies and wind 
opening angles that varied with SMBH accretion rate 
such that both quantities are small at small accretion 
rate and reach a specified maximum at the Eddington 
rate. 

Independently of the mechanical feedback model, the 
radiative luminosity of the AGN is given by 



L = cem-Mbhc 



(18) 



where the electromagnetic efficiency is gi ven by the ad- 
vection dominated accretion flow inspired (N aravan fc Yil 
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U99l formula: 



e EM — 



eoAm 
1 + Am ' 



(19) 



and A = 100 and eo = 0.1. The dimensionless mass 
accretion rate is 



Af BH _ e M BH c 2 
M ¥jM ^Edd 



(20) 



where LEdd is the Eddington luminosity. 

The mechanical wind efficiency is given by ew = con- 
stant for A models. The specific parameters for the sim- 
ulations presented here are given in Table [TJ The present 
work considers two-dimensional analogs of the simpler A 
models only. 

We also ran tw o-dimensional analogues of the 
ICiotti et al.l (|2009bf) B models where the mechanical 
feedback efficiency falls off at low Eddington ratios 
(as expected from computational studies of the cre- 
ation of BLWs in the inner few hundred gravitational 
radii (Kurosawa fc Progal l2009at Kurosawa et al.l [20091 : 
Kurosawa fc P roga 2009b)). In one dimension, mechani- 
cal feedback in the B models was sufficient to regulate the 
SMBH growth. However, in two dimensions equivalent 
models were found to undergo too much SMBH growth 
and had distributions of Eddington ratios very differ- 
ent from those observed for actual galaxies. Comparing 
one-dimensional and two-dimensional A models reveals 
that feedback is in less effective in two dimensions: more 
energy is required for the SMBH to effectively regulate 
its growth. For two dimensions, the low efficiencies at 
low Eddington ratio combined with the diminished ef- 
fectiveness of feedback in two dimensions make the two- 
dimensional B models observationally unacceptable. 

The model most similar to that used by 
iDi Matteo et al.l (|200ai is A0. The model has an 
efficiency of 5 x 10~ 3 , but we take into account the 
fact that the wind carries not only energy but also 
momentum and mass into the ISM. Matching our 
mome ntum input rate to that used by iDebuhr et al.l 
(2011) gives e w = v w € t/2c ~ 0.04, where r = 25 
is their chosen wavelength-averaged optical depth 
parameter and v w = 10, 000 km s _1 is our chosen BAL 
wind velocity parameter. The present models do not 
include a simulation with such a high feedback efficiency. 
Our model that m ost closely approximates those of 
IDebuhr et al.l ([20111 ) is again A0, with the caveat that 
we inject l ess momentum for th e same BH luminosity. 
Note that [Debuhr et al.l ([20111 ) effectively assume a 
much lower wind velocity (not intended to model BAL 
winds) and as a result the energy actually injected into 
th e ISM in their s i mulat ions is similar to that used 
in IDi Matteo et al.l (|2005| ). The above considerations 
are intended to determine what parameters we would 
need to adopt for our present simulations in order to 
produce models similar to existing ones by other authors 
given our assumed BAL velocity. A closer study of the 
similarities and differences between these AGN feedback 
models is the subject of a forthcoming paper. 

2.3. Angular Momentum 

The one-dimensional simulations did not permit the 
simulated gas to have non-zero angular momentum. The 



two-dimensional simulations assume axisymmetry and 
compute the velocity in the 0-direction. We must assume 
an angular momentum profile for the injected gas. In 
the present simulations, we avoid forming a rotationally 
supported gas disk by choosing the radius of centrifugal 
support to be inside the innermost grid cell. 

This assumption gives rotation velocities compa- 
rable to the slowest of th e SAURON slow-rotators 
(|Emsellem et al.l 12004 120071 ). Thus, our model galaxy 
is representative of a rather small set of observed galax- 
ies. However, the advantage of assuming very little 
rotation is that it allows us to avoid specifying an 
ad hoc prescription for angular momentum transport. 
The low angular momentum case also allows the clos- 
est comparison to the existing one-dimensional simula- 
tions, allowing us to isolate the effect of dimensionality. 
In future work we plan to implement standard recipes 
for angular momentum tr ansport, including an a disk 
(|Shakura fc Sunvaevl 119731) and models based on grav- 
itational torques ([Hopkins fc Quataertl 12010a!) . Taking 
angular momentum transport into account is likely to 
have a significant effect on our computed results. Our 
consideration of the low angular momentum case will fa- 
cilitate isolating the effect of angular momentum in fu- 
ture analysis. 

The assumed net specific angular momentum of the 
stars providing gas in the simulation is thus taken to be 
given by 

1 d 1 R 

v^(R) oR fa j 

where v<p is the velocity in the azimuthal direction, 
R is the distance to the z-axis, a is the central 
one-dimensional line-of-sight velocity dispersion for the 
galaxy model, and d, /, and j are adjustable parame- 
ters controlling the angular momentum profile at small, 
intermediate, and large radii, respectively. This param- 
eterization gives solid body rotation at radii less than d 
and constant specific angular momentum of j at large 
radii. The second term on the right prevents the rota- 
tional velocity from exceeding fa at any radius. If there 
is a range of values of R for which the fa term dominates, 
then is a region of intermediate radii with constant rota- 
tional velocity. 

2.4. Initial Conditions 

The initial conditions are chosen to match those 
of the one-dimens ional simulations, fully described in 
ICiotti et al.l ([2010D . Briefly, the total gravitational po- 
tential is that of a singular isothermal sphere with one- 
dimensional velocity dispersion 260 km s _1 and with a 
central point mass with initial mass 3 x 10 8 M Q . The 
stellar distribution is a Jaffe profile with projected half- 
light radius 6.9 kpc, mass-to-light ratio 5.8 in solar 
units, and total mass 3 x 10 11 M Q . These parame- 
ters are chose n to put the initial galaxy on the Funda- 
ment al Plane ([Diorgovski fc pavis]ll987t iDressler et al.l 
19871). the Faber-Jackson (119761 ) relation, and the 



Magorrian et al.l (|1998t ) relation. All of the relevant dy- 



namical properties of the galaxy models are given in 
ICiotti et al.l J2009a) . The gas density is initially set to a 
very low value so that the gas in the simulation comes al- 
most exclusively from explicit source terms arising from 
stellar evolution. 
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2.5. Boundary Conditions 

We assume reflecting boundary conditions on the 6 
boundaries occurring at either pole. On the outer radial 
boundary we assume that all fluid quantities are con- 
stant. This allows both outflow and inflow depending on 
the state of gas just inside the outer boundary. 

On the inner radial boundary we assume reflecting 
boundary conditions if the innermost radial velocity is 
positive. Physical processes that inject mass (such as 
the BAL wind) are handled as explicit source terms act- 
ing in the first set of radial cells. This allows us to easily 
handle the cases of very strong, very weak, and inter- 
mediate flows of energy, mass, and momentum onto the 
computational grid with a single source term. If the in- 
nermost radial velocity is negative, we use an outflow (off 
of the computational grid toward the center of the sim- 
ulation) boundary condition where all fluid variables are 
constant across the boundary. The one exception in this 
case is the radial velocity itself. If the cell is inside the 
locally estimated Bondi radius (that is, the Bondi radius 
has been resolved), then no limit is imposed upon the 
inflow velocity. However, if the cell is outside the locally 
estimated Bondi radius (the Bondi radius is unresolved) , 
then the radial velocity is limited to the velocity that 
would result in mass transport consistent with the Bondi 
accretion rate. That is 



v r ,o 




if v r ,x > 

if i> r> i < 0, r x < r B 

if V r ,x < 0, r x > r B 



(22) 



As discussed above, we are careful to resolve the Bondi 
radius even for gas heated to the maximum expected 
temperature for infalling gas, the Compton temperature. 
However, the BAL wind has an even higher specific en- 
ergy. Typically, the BAL wind is flowing strongly out- 
ward and the Bondi radius is always resolved. However, 
it is possible for the wind to fill the inner region of the 
galaxy with gas heated above the Compton temperature. 
In this case, Equation (j2"2")l ensures that the SMBH ac- 
cretion rate is reduced accordingly. 

It is unfortunately not possible to cleanly separate the 
SMBH scales from the galactic ones. Physically, this 
is because the source of gas in the galaxy is the stars, 
evolving on galactic timescales. Meanwhile, the SMBH is 
easily able to affect gas on kiloparsec scales. This ties the 
scales together in a feedback loop that cannot be easily 
separated into two separate simulations, or a simulation 
plus a sub-grid model. Furthermore, the classic Bondi 
solution is for the simple case of a point mass. If the 
mass enclosed by the Bondi radius is dominated by the 
galaxy rather than the central BH, then the classic Bondi 
solution is modified, particularly if the galactic mass is 
not spherically symmetric. 

3. RESULTS 

The character of the SMBH accretion changes dramat- 
ically in going from one to two dimensions. In one di- 
mension, accretion events happen when a cold shell of 
gas forms at ~ 100 pc. The shell falls into the SMBH 
as a unit and, after a series of sub-bursts in which di- 
rect and reflected shock waves interact and carry new 
material for accretion on the SMBH, triggers a dramatic 
release of energy from the SMBH. This leaves a sphere 



of hot gas at the center of the simulated galaxy. Subse- 
quent cold shells can only reach the center when the gas 
beneath them either cools or is compressed inside the in- 
nermost grid point. The hot gas generated by radiative 
and mechanical AGN feedback is able to prevent SMBH 
accretion until it cools. These processes lead to dramatic 
bursts followed by long, extremely quiet periods, spaced 
by the cooling time of the central gas. As the collapsing 
cold shell sits on top of hot and low density gas, it was 
already clear from the previous one-dimensional simula- 
tions that the shells are RT unstable. 

In two dimensions, cold gas forms again at ~ 100 pc. 
However, cold gas takes the form of rings rather than 
shells due to the classical RT instability. Shells some- 
times form, but they fragment quickly. If there is hot gas 
beneath the cold ring, both the RT and convective insta- 
bilities operate to allow the cold gas to fall to the center 
of the simulation unimpeded. These instabilities cannot 
operate in one dimension. In two dimensions, they allow 
both higher and more chaotic SMBH accretion rates. 

Figure [T] shows the simulation during a relatively quiet 
period. Figure [2] shows the start of an accretion event. 
A cold blob of gas is freely falling to the center of a 
two-dimensional simulation. Figure [3] shows a simulation 
snapshot just after an accretion event with the bipolar 
BAL wind flowing away from the center of the simula- 
tion. Gas is able to continue to fall into the SMBH via 
the simulation midplane. Finally, Figure |4] shows a simu- 
lation snapshot significantly after an accretion event (but 
not so long that the galaxy is able to return to a quiescent 
state). Dense overlying gas has caused the BAL wind to 
become nearly isotropic, making additional accretion via 
the midplane impossible. 

3.1. Black Hole Growth 

Figure [5] shows SMBH mass versus mechanical feed- 
back efficiency for one- and two-dimensional A models 
as well as one-dimensional B models. The SMBHs un- 
dergo more growth in two dimensions at a fixed feedback 
efficiency owing to RT instabilities, making it easier for 
cold gas to fall in. 

3.2. Time Dependence of SMBH Accretion 

The character of the time dependence of SMBH accre- 
tion is much different in two dimensions than in the one- 
dimensional case. In one dimension, SMBH accretion 
occurs in a few bursts well separated in time. Occasion- 
ally a given burst will have a complex character, being 
composed of many sub-bursts. In two dimensions, there 
are still occasionally events that can be characterized as 
bursts followed by quiescent periods. However, the qui- 
escent periods are shorter and the SMBH is more active 
during them compared to the one-dimensional case. Fur- 
thermore, during times far from a major burst when gas 
is building up in the galaxy, the SMBH is far more active 
than in the one-dimensional case and the accretion takes 
on a stochastic character. 

Figure [6] shows the Eddington ratio L/L^dd as a func- 
tion of time for part of the A2 simulation. The SMBH ac- 
cretion is much more chaotic than in the one-dimensional 
case. In the one-dimensional case, a cold shell forms and 
falls in as a unit only after the gas interior to it has also 
cooled. In two dimensions, these same cold shells form, 
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Fig. 1. — Simulation snapshot during a quiescent period. On the left, gas density in number of protons per cubic centimeter. In the center, 
log sound speed in kilometers per second. On the right, the radial velocity in kilometers per second. The x- and y-axes are logarithmic 
in the distance to the SMBH. Outflowing gas from past accretion events is visible from 300 pc to 10 kpc. Persistent low-level accretion 
maintains a narrow continuous outflow near the poles. 
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Fig. 2. — Snapshot from an axisymmetric simulation showing a cold blob falling to the center of the galaxy. The cold gas was produced by 
enhanced cooling in an overdense quasi-spherical shell with a covering fraction of about one-third of the sphere. The gas quickly collapses 
to a ring with a small covering fraction and/or fragments as it freely falls to the center of the simulation. 
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Fig. 3. — Simulation snapshot during an accretion event. A significant quantity of hot, outflowing gas injects energy and momentum into 
the interstellar medium at r ~ 1 kpc. 
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Fig. 4. — Simulation snapshot showing the final stages of a major accretion event. A hot, expanding bubble of gas extends to 100 pc, 
shutting down further SMBH accretion. Dense, overlying gas has caused the initially bipolar BAL wind to become quasi-spherical. The 
hot bubble is breaking through the overlying gas in at the north pole, which will lead to a unipolar wind. 
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Fig. 5. — Final SMBH masses vs. mass- averaged wind Efficiency for one- and two-dimensional A models. The low-efficiency two- 
dimensional models produce SMBHs at the upper end of the range of observed central SMBHs for the characteristics of our model galaxy. 
Two-dimensional models allow more SMBH growth at a given efficiency because instabilities allow cold shells of gas to break up and fall 
into the center of the simulation with greater ease than in the one-dimensional case. A more accurate two-dimensional treatment of the 
radiative transfer in optically thick regions would probably lead to a reduction in the SMBH masses. 

but they immediately break up into blobs that are RT 
unstable. The cold blobs fall in to the SMBH on a free- 
fall time in a much more disorganized fashion compared 
to one dimension. 

Figure [7] shows the power spectrum of the Eddington 
ratio as a function of time for the A2 simulation. In spite 
of the complex, detailed implementation of the physics of 
stellar evolution, atomic cooling, star formation, super- 
novae, and AGN feedback, the power spectrum reveals 
a simple, gently sloped power law (roughly frequency to 
the one-fourth power) at low frequencies. At high fre- 
quencies, the power spectrum falls off as 1// because of 
the filtering effect of the SMBH accretion disk. 

Figure [8] shows the cumulative distribution of Edding- 
ton ratios for the A2 simulation. The black (red) lines 
show the cumulative time the simulation spends above 
(below) the given Eddington ratio. Solid lines show the 
cumulative distribution of Eddington ratios for the entire 
simulation, while the dashed lines show the distribution 
for the final 2 Gyr (corresponding to low redshifts). Low- 
redshift observational constraints for the fraction of 3 x 
1O 8 M BHs accret ing at 1% and 10% of t he Eddington 
rate a re tak e n from iHeckman et al.l (l2004D.lGreene &; Hoi 
(12001 . iHol d200l . and TKauffmann fc Heckmanl (12009^ 
A constraint on the fraction of high-redshift Lyman- 
break galaxi e s sho wing nuclear activity is taken from 
ISteidelet all (p00l . 

There is reasonable agreement between the simulations 
and observational constraints for a mechanical efficien- 
cies €w between 10 -3 and IP" - 1 . This range is in agree- 
ment with the recent study by lArav et alJ (|2011f >. They 
presented new observations as well as values drawn from 
the literature for E and M for five broad-line systems. 
The mean value of loge^ for the five systems is —4.2 
with a mean error of 0.2 dex for each system and a stan- 
dard deviation of 0.4 dex for the five points. 

3.3. Disk Wind Opening Angle 



The opening angle of the disk wind is expected to have 
an effect on the SMBH growth by controlling the effec- 
tiveness of the AGN feedback. At very small opening 
angles, the energy injected by the SMBH drills through 
the nearby gas and exits the central region as a narrow 
beam. All of the energy is eventually deposited at radii 
beyond one kpc, leaving the cooling and infalling gas in 
the central region essentially undisturbed. This should 
result in large SMBH growth. 

As the opening angle increases, one might expect that 
mechanical feedback would become more and more ef- 
fective as the wind interacts with an increasingly large 
fraction of the gas surrounding the SMBH. This would 
imply that isotropic feedback would result in the least 
SMBH growth. 

Figure [5] shows final SMBH masses versus wind open- 
ing angle. The largest SMBH growth occurs when the 
wind is isotropic. A very wide opening angle allows the 
wind to effectively couple to all of the gas at small ra- 
dius, temporarily suppressing SMBH accretion. How- 
ever, the wind's outward progress is soon stopped be- 
cause the wind must essentially lift the entire overlying 
gaseous atmosphere. The wind is not able to move a 
significant amount of gas to a sufficiently large radius to 
prevent eventual accretion onto the SMBH. 

This implies that there is an opening angle that re- 
sults in minimal SMBH growth. The angle will likely be 
a function of the mechanical feedback efficiency as well 
as the other physical parameters for the gravitational po- 
tential, stellar distribution, mass-loss rates, and so forth. 
Figure |9] shows the increase in SMBH growth at large 
opening angles, but does not show evidence for the ex- 
pected increase in SMBH growth as the opening angle 
becomes very small. However, it is worth noting that 
the change in BH mass growth is less than a factor of 
two for a wide range of wind opening angles. 
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Fig. 6. — Eddington ratio as a function of time, for three different time intervals in the A2 simulation. 



3.4. Star Formation and Galactic Winds, and Gas 
Content 

Table [1] gives total mass of stars formed, total mass 
of gas driven beyond 10i? e , and final mass of gas 
within 10R e for each simulation. Star formation con- 
sumes about 30% of the total mass budget in the two- 
dimensional simulations and is very insensitive to the 
details of the AGN feedback. This is in good agreement 
with the one-dimensional simulations at low mechanical 
efficiencies. However, at high feedback efficiencies, the 
one-dimensional simulations drive significant quantities 
of gas out of the galaxy, leading to low star formation 
rates and low final gas content. In this respect the one- 
dimensional and two-dimensional simulations disagree. 
However, this is to be expected since assuming spherical 
symmetry gives the most favorable situation for turning 
a central energy source into a global outflow. In two 
dimensions, energy can escape via low-density channels 
and fail to participate in driving an outflow. 



Figure [10] shows the mean mechanical energy input 
versus the mean efficiency for one-dimensional and two- 
dimensional A models. For two-dimensional A models, 
the energy input is nearly constant — the SMBH accre- 
tion self-regulates to provide energy at this rate. The 
one-dimensional A models have lower energy input rates. 
That is, two-dimensional models require more energy to 
reach equilibrium between inflow (due to cooling) and 
outflow (due to mechanical feedback). 

4. CONCLUSIONS 

We have performed two-dimensional simulations of the 
entire cosmic history (12 Gyr) of an isolated L* elliptical 
galaxy. Planetary nebulae and red giant winds produced 
by evolving low-mass stars serve as the source of gas in 
the galaxy. This gas finally ends up either in the central 
BH, in long-lived low-mass stars (formed in the simula- 
tion), in the ISM within the galaxy (at the end of the 
simulation), or outside the galaxy as part of the inter- 
galactic medium. As gas finds its way to one of those four 
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Fig. 7. — Power spectrum of the dimensionless mass accretion rate M-gii/M^d f° r A2. The units of the y-axis are arbitrary. For 
frequencies higher than the inverse of several times the central accretion disk timescale, the power spectrum is proportional to 1//, the 
frequency response of the low-pass filter applied by our central accretion disk. For lower frequencies, the power spectra are determined 
by the physics of gas input, cooling, and feedback. At these low frequencies, the power spectrum is nearly flat: oc / _1 ' 4 . That is, the 
accretion onto the SMBH has a power spectrum nearly the same as white noise. 



TABLE 1 
Properties of Computed Models 



Model 


c\v 


(cem) 


log AA/bh 


log 


AM, 


log AM W 


log 


J "gas 


AO 


5 x icr 3 


0.041 


7.40 




9.98 


10.28 




9.67 


A05 


icr 3 


0.060 


8.10 




9.94 


10.33 




9.50 


Al 


2.5 x 10" 4 


0.067 


8.71 




9.88 


10.35 




9.45 


A2 


1(T 4 


0.067 


9.18 




9.91 


10.29 




9.49 


A3 


5 x icr 5 


0.063 


9.55 




9.94 


10.17 




9.58 


A4 


3 x icr 5 


0.060 


9.72 




9.93 


10.11 




9.50 


A5 


l x icr 5 


0.048 


9.71 




9.98 


10.03 




9.50 


A6 


3 x icr 6 


0.047 


9.69 




10.00 


10.03 




9.52 


A7 


l x icr 6 


0.047 


9.69 




9.99 


10.03 




9.52 



Note. — Final properties of the simulated galaxies, where eyy is the me- 
chanical wind efficiency, (cbm) is the mass- weighted mean radiative efficiency, 
AAfBH i s the change in the mass of the BH, AM* is the mass of stars formed 
during the simulation, AMy is the total gas mass driven beyond 10 effective 
radii, and logA/ gas is the total mass remaining within 10 effective radii. 



final states, it can engage in star formation, mass/energy 
injection into the ISM via Type la or Type II supernovae, 
or radiative cooling. 

The primary purpose of the code is to implement the 
well-developed physic al model of mechanic al and radia- 
tive AGN feedback of iCiotti et ail (j2009bD . This model 
includes momentum and energy imparted to the gas by 
Compton scattering and atomic lines via radiation from 
the central SMBH. It also includes a self-consistent model 
of mechanical feedback via a broad-line wind. The pri- 
mary difference between the two-dimensional code and 
the one-dimensional code (upon which it is based) is that 
the two-dimensional code does not yet include the physics 
of optically thick radiation pressure on dust grains. Ex- 
periments with the one-dimensional code have shown 
that turning off the dust opacity does not greatly affect 
the SMBH growth, although other aspects of the simula- 
tion results such as the star formation rate are affected. 



Multi-dimensional simulations are a significant im- 
provement over the one-dimensional simulations because 
the additional degrees of freedom allow classical instabil- 
ities such as the RT or Kelvin-Helmholtz instabilities to 
operate. The RT instability in particular plays an im- 
portant role because AGN feedback tends to produce a 
hot bubble of gas near the SMBH, effectively halting ac- 
cretion. In one dimension, the bubble cannot move away 
from the SMBH and continues to suppress subsequent 
infall until it itself radiatively cools. In two dimensions, 
the hot bubble buoyantly moves up through the ISM and 
additional cold gas generated at radii of 0.1-1 kpc can 
freely fall into the center of the simulation. 

In this work, we have assumed a low specific angular 
momentum profile for the model galaxy in order to match 
the one-dimensional simulations as closely as possible. 
The effect of angular momentum transport in galaxies 
with more typical specific angular momenta is likely to 
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Fig. 8. — AGN duty cycle as a function of Eddington ratio for the A2 simulation. The black/red lines show the cumulative time 
above/below the given Eddington ratio. Solid lines are computed using the entire simulation time, while dashed lines use only the final 
2 Gyr. Points are observational constraints and are similarly colored black or red according to wh ether they are measur ements of the 
fraction of ob j ects a bove or below the given Eddinto n ratio. Downward-poin t ing tr iangles are from Hcckman ct al. (2004), circles from 
I Greene fc Ho (2007), upward-pointing triangles from Kauffmann & Hcckman (2009), and squares are from Ho (2009). All of these use 
low-redshift observations and should thus be compared to the dashed lines. The star is a constraint from Stcidcl ct al. (2003) for the 
Lyman-break galaxies at high redshift showing nuclear activity, and therefore should be compared to the solid lines. There is significant 
observational disagreement about the fraction of 3 X 10 8 Mq BHs accreting at 1% or 10% of the Eddington rate, even at low redshift. 
Nevertheless, the A05 and A2 simulations fit both the high-redshift and low-redshift constraints reasonably well given the state of the 
observational data. 

be very important. We plan to implement angular mo- 
mentum transport in future work in order to isolate its 
affect on our results. 

An important feature of our simulations is that we have 
taken care to resolve the Bondi radius even for gas heated 
to the Compton temperature by the AGN. The gas heat- 
ing terms due to AGN radiation go as l/r 2 , becoming 
much more effective as a given parcel of gas moves to- 
ward the BH. If the Bondi radius is resolved even for the 
hot gas, then our conclusion about whether or not a given 
parcel of gas makes it to the center of the simulation is 
secure. If the Bondi radius is not resolved, then radia- 
tive heating that would have occurred between the inner 
radius of the simulation and the physical Bondi radius 
might have increased the temperature of the gas to the 
point that thermal energy dominated over gravitational 
energy — the gas would not make it to the SMBH. Thus 
if the Bondi radius is not resolved for the hot as well as 
the cold gas, we cannot be sure whether or not a given 
parcel of gas actually makes it to the SMBH. 

Resolving the Bondi radius for the hot gas requires high 



central spatial resolution. The large velocities of the BAL 
wind entering the grid near the center necessitate time 
steps of order one year. The timescale relevant for stellar 
evolution, the physical source of the gas, is several Gyr. 
We have opted to use rather coarse spatial resolution in 
order to run the simulations long enough to see evolution 
over the entire cosmic history of the galaxy in question. 
The primary differences between the previously pre- 
sented one-dimensional simulations and the present two- 
dimensional simulations are that the character of the ac- 
cretion as a function of time changes from well-separated, 
dramatic bursts (one dimension) to chaotic, stochastic 
accretion (two dimensions). There are two reasons for 
this. First, the cold, infalling gas generated at a few hun- 
dred parsecs tends to fragment and fall into the SMBH as 
several discrete blobs (two dimensions) rather than as a 
coherent spherical shell (one dimension). Second, the RT 
instability allows hot gas generated via radiative heating 
by the AGN to move out of the way, allowing subsequent 
blobs of cold gas to fall into the SMBH unimpeded. Thus, 
a burst of accretion in two dimensions is much less effec- 
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Fig. 9. — Change in SMBH mass vs. wind opening angle for eyy = 10 — 3 . Mechanical feedback becomes less effective for opening angles 
larger than 45°, corresponding to a covering fraction of 1/4 of the sky. Overall, the change in BH growth is less than a factor of two for a 
wide range of wind opening angles. 
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Fig. 10. — Time-averaged mean energy input (f Edt/ j dt) via AGN mechanical feedback vs. mass-averaged mean feedback efficiency 
(f eradt/ J rhdt). For the two-dimensional simulations, at the low-efficiency end, energy input falls off because the SMBH is consuming 
all of the mass available to it: AMbh is at a maximum but ew is still falling. At the high-efficiency end, energy input approaches a 
constant — the simulation self-regulates to an energy input of ~ 5 X 10 41 erg s — 1 . This implies a nearly perfect, inverse relationship between 
AMbh an< f e W a *- i- ne high-efficiency end. The one-dimensional simulations have smaller energy inputs at a given efficiency, reflecting their 
smaller SMBH growth rates. 



five at stopping subsequent accretion than is the case in 
one dimension. 

Another way of saying this is that mechanical feedback 
is less effective in two dimensions than one dimension: 
it takes more energy for a given galaxy model to reach 
equilibrium between cooling-driven inflow and feedback- 
driven outflow/heating. At a fixed mechanical efficiency, 
the SMBH self-regulates to larger mean accretion rates 
in order to provide the required additional energy. 

Two-dimensional outflows are less effective than spher- 
ical ones in two critical aspects, less effective in protect- 
ing the SMBH from infalling gas and less effective in 



ejecting gas from the galaxy. 

For our models with wind efficiencies of 10 -3 and 
10 -4 , the distribution of Eddington ratios is in reason- 
able agreement with the observed fraction of galaxies 
above/below a given Eddington ratio over four orders of 
magnitude in luminosity. These values of the efficiency 
are both phy sically plausible and supported by recent 
observations (|Arav et al.ll201ll ). 

Our plans for future work include implementing opti- 
cally thick radiative transfer on dust to bring the physics 
of our treatment of AGN and star-formation feedback up 
to the state of the art. Furthermore, the B models for 
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mechanical feedback are more realistic in the sense that 
the efficiency and opening angle both become small at 
low accretion rates, in line with observational and the- 
oretical expectations. The one-dimensional simulations 
produced B models in good agreement with observations 
(jCiotti et al.l 120 10). but using identical parameters for 
two-dimensional simulations did not result in the same 
good agreement. A parameter study of B models in two 
dimensions is necessary to identity viable ranges of the 
parameters. Finally, there remain many potentially im- 
portant physical processes that we have not attempted to 
model in this work — energy transport via conduction and 
the extent to which conduction is suppressed by tangled 
magnetic fields among them. 

In this paper, we have focused on the comparison be- 
tween one-dimensional and two-dimensional AGN feed- 
back models with nearly the same implemented physics 
and spatial resolution. A future study will put this 



model into better context by comparing more thor- 
oughly to t he existing AGN feedback models in the liter- 
ature (e.g.. iDi Matteo et alj|2005t Uohansson et al.H2009t 
IDebuhr et al.ll2010l h 
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